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Eikonal, or ray tracing, methods are commonly used to estimate the propagation of radio frequency fields 
in plasmas. While the information gained from the rays is quite useful, an approximate solution for the fields 
would also be valuable, e.g., for comparison to full wave simulations. Such approximations are often difficult 
to perform numerically, because of the special care which must be taken to correctly reconstruct the fields near 
reflection and focusing caustics. In this paper, we compare the standard eikonal method for approximating fields 
to a method based on the dynamics of wave packets. We compare the approximations resulting from these two 
methods to the analytical solution for a lower hybrid wave reflecting from a cutoff. The algorithm based on 
wave packets has the advantage that it can correctly deal with caustics, without any special treatment. 



I. INTRODUCTION 

Eikonal methods — also known as ray tracing or Wentzel- 
Kramers-Brillouin (WKB) methods — are often used to model 
the propagation and absorption of radio frequency waves in 
plasmas (for ray tracing applied to lower hybrid waves see O)- 
ISQ). The results obtained from ray tracing can be compared to 
full wave simulations, in order to reproduce, supplement, and 
better understand those simulations |@]. Since the ray trac- 
ing method was devised as an approximate method for prob- 
lems where the wavelength is short compared to gradients in 
the medium, there are problems which cannot be accurately 
solved using ray tracing. There are also cases where diffrac- 
tive effects like focusing and reflections from cutoffs are not 
handled well by traditional ray tracing methods. These caus- 
tics are often interpreted as a breakdown of the eikonal ap- 
proximation used to justify ray tracing, and therefore an in- 
dication that the ray tracing results are not valid f?]. How- 
ever, more advanced eikonal methods have been developed 
which allow ray tracing to be used even in the presence of 
caustics. These methods, based primarily on the work of 
Maslov 1 8], are routinely and successfully applied to wave 
problems in various fields, including atomic, molecular, and 
optical (AMO) physics |9, 10] and geophysics |IIl[ll]. Be- 
cause these advanced eikonal methods are able to deal with 
caustics, it becomes possible to construct approximate wave 
fields based on ray tracing data. 

However, the techniques based on the work of Maslov re- 
quire careful application, and special attention must be paid to 
the regions around caustics. Different types of caustics require 
various methods to deal with them, and it becomes a "labor of 
love"[i l3J to perform the field reconstruction. Because of the 
difficulty of properly applying the Maslov methods, it is quite 
challenging to create an numerical algorithm for field recon- 
struction based on these techniques. 

In this paper we have, for the first time that we are aware 
of, employed these methods to reconstruct RF wave fields 
in situations relevant to realistic fusion plasmas. The prob- 
lem of a lower hybrid wave reflecting from the density cut- 
off is considered. Given certain approximations, this problem 
can be solved analytically for a cold plasma slab model. We 
first describe the model system and the analytical solutions, 
which can be compared to the approximate solutions that will 
be presented. We then use the standard Maslov technique to 



construct an eikonal solution which is valid even in the re- 
gion near the caustic. This will illustrate the difficulties that 
arise when trying to use standard eikonal methods to automat- 
ically construct approximate solutions numerically. We then 
describe a technique for constructing fields that is based on an 
approximation for wave packet dynamics. Application of this 
method to the lower hybrid cutoff problem shows how it is 
well suited for use as a numerical algorithm for the construc- 
tion of approximate solutions, even in the presence of caus- 
tics. The paper ends with a comparison of the two approxi- 
mate solutions to the analytical solution, and a discussion of 
the results. 



II. THE LOWER HYBRID CUTOFF 

In this section, we will describe the problem of a lower hy- 
brid wave reflecting from the density cutoff, and derive the 
analytical solution. In Sec. |Vl we will use this solution as a 
benchmark with which to compare our eikonal solutions. 

As in Stix llT4fl . the equation we want to solve is 

-Vx (VxE) + ^K-E = 0, (1) 

which is the wave equation for the electric field E in a plasma 
with dielectric tensor K. The time dependence in this equation 
has already been removed through the substitution idt — > CO. 
This wave equation can be written in a compact form as 

D(x,V).E(x)=0, (2) 

where D is the wave operator, which is a function of the non- 
commuting operators x and V. The ambiguity in the ordering 
of these operators is resolved by using the Wigner-Weyl for- 
malism, a description of which can be found in Appendix B 
of Ref. jl5i]. This formalism provides a well defined process 
for computing the dispersion matrix D(x,k) associated with 
an arbitrary wave operator D. However, for the problem to 
be considered here, this process reduces to making the substi- 
tution —/V — > k = Nco/c, since the ordering of the multipli- 
cation and gradient operators is unambiguous for this model. 
We will now give the detailed description of our model and its 
dispersion matrix. 
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A. Description of the model 

Consider a two-dimensional (2D) cold plasma slab model 
with a linear density gradient n{x,z) = n'x. The magnetic 
field is uniform, and points in the z direction: B(x,z) ~ Bqz. 
In this paper, we will use the parameters Bq = 5.5 T, and 
n' — 3 X lO'^m^"*, / — 4.6GHz. We will use the cold plasma 
dispersion matrix from Stix 1141 : 
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and cOpj and Clj are, respectively, the plasma and cyclotron 
frequencies of the particle species. Also, Ej is the sign of 
the charge of the species. The frequencies are given by 
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Since in this model the magnetic field is assumed to be con- 
stant in space while the density is allowed to vary, all of the 
spatial variation in D is contained in the plasma frequencies 

COpj. 

In this paper, we are interested in the reflection of a wave 
from the slow wave cutoff at P{x) = 0. In this region, the 
spatial variation which we are most interested in is the vari- 
ation in P{x). Here, S{x) ~ 1 and D{x) ~ 0. In order to be 
able to derive an analytical solution to our problem, we will 
modify the cold plasma dispersion relation by setting D = Q 
and 5=1. While this does change the problem, it retains the 
physics of interest [the variation in P{x) which gives the slow 
wave cutoff], and it allows us to derive an exact solution that 
we can compare to our numerical field reconstructions given 
in the later sections. These modifications (D — and 5=1) 
will be retained throughout the remainder of the paper. 

With these modifications, the dispersion function for our 
model becomes: 
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FIG. 1: Dispersion function for the slab model used in this paper; 

Af,, = 0, A', = 2. 



B. The analytical solution 

In this paper we studying the reflection of a lower hybrid 
wave from the slow wave cutoff. For such a reflection, the 
dominant spatial variation is in the direction perpendicular to 
the cutoff (the x direction), so we have modeled our plasma as 
being uniform in the z direction. Since there is no variation in 
the z direction, we can Fourier transform Eqn. ([T]) in z, giving 
us a set of equations for 'E,{x;Nz). These equations can also be 
obtained from the dispersion matrix in Eqn. ^ by making the 
substitution Nx = ^i{c/ (o)dx, which gives 
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Once these equations are solved for each A'^, we can form a 
general solution to the problem by summing each A^, mode 
with some arbitrary spectral density F{Ny): 



E(x,z) 
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We can now solve for each Fourier mode independently, and 
we obtain a solution which can be written in terms of the Airy 
function Ai and its derivative Ai'. The vector components of 
the solution are: 
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Solving for S^{x,z,Nx,N^) = gives the dispersion relation for 
this problem. Figure [T] shows the dispersion function for the 

case Ny ^0,N.^ 2. 
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FIG. 2: Spectrum in A^^ for the analytical solution. 
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FIG. 3: (Color online) Absolute value of for the analytical solu- 
tion. The black line is the line of stationary phase given in Eqn. dlVb . 
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Note that for our linear density profile, {d^P) is a negative 
constant, and therefore a is a constant. Also, since we will 
take |A^^| > 1, we have that a is a positive real number In 
order to create a beam-like solution, we choose F{N^) to be a 

gaussian, with its peak at Ni^\ and a width of Gn.'- 
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We can get an idea of what the solution with this spectral 
density will look like by examining the integral in Eqn. ( fT2b . 
Since the spectrum is peaked (see Fig. |2]i, we can develop 
some intuition about the solution by applying the stationary 
phase approximation to the integral. Rather than using the 
stationary phase method to obtain an analytical approximation 
to the solution, we will simply use it to find the curve in the 
(ji:,z)-plane where we expect the maximum field amplitude. 
This curve will describe the center of the resulting "beam." 
Using the asymptotic expressions for the Airy functions lfT6ll . 
we obtain the equation for the curve: 
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This curve, along with the amplitude of the £- component of 
the solution, are shown in Fig.|3] 



III. STANDARD EIKONAL FIELD RECONSTRUCTION 

In this section, we apply the standard eikonal approxima- 
tion to the problem of the reflection of a lower hybrid wave 



from the cutoff. We use the method detailed by Maslov and 
Fedoriuk in Ref . fs") . Rather than giving detailed derivations 
of the equations used for this approximation, we will simply 
give the relevant equations and explain how they are used. For 
more details, the interested reader is encouraged to refer to 
Ref. js] or Ref. [10], where this method is applied to a quan- 
tum scattering problem. 

In the standard eikonal method, a family of rays is used 
to construct a field, which is an approximate solution to the 
wave equation being studied. This approximate solution has 
the eikonal form 



E(x,f) =A(x)e(x)e'®W-''"'. 



(18) 



The family of rays traces out a surface in (x,k) phase space, 
called a Lagrange manifold. The family of parameterized rays 
forms a coordinate system on this surface, with the ray "time" 
T as one coordinate direction, and the ray label j3 forming the 
other coordinate. The individual rays are found by solving 
Hamilton's equations for the ray position 



H{x,k), 



dx 

dx dk 

dk d 

— = H[x,k), 

dx dx ^ ' 



(19) 
(20) 



where H{x,k) is a zero eigenvalue of the dispersion matrix 
D(x,k), with eigenvector e(x,k). 

The spatial gradient of the phase is given by V0(x) = k, 
where x and k are evaluated along the ray. Given an initial 
phase 0(xo) = ©o defined by boundary conditions at xq we 
can integrate this equation to find 0(x): 



0(x) = 00+ / k(xVx' 



(21) 
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Here, k(x) refers to the function found by inverting the ray 
data (x(T,/3),k(T,j3)) to obtain k as a function of x along the 
ray. In practice, the ray parameter T and label j3 are often used 
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for integrating the equation for the phase, since that is simpler 
than solving for k(x). In this case, we have 

, . , , rPf 
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It turns out that — for single-valued k(x) — the result of the in- 
tegral in Eqn. ( |2TI ) is independent of the path of integration. 
Because of this, the limits on the integrals in Eqn. ( l22b need 
only be chosen so that x(t,-,/3,) = xq and x(t/, j3/) is the point 
where 0(x) is to be evaluated. 

This integration must be treated slightly more carefully 
when k(x) is double (or multiply) valued. This situation oc- 
curs when Lagrange manifold folds back over itself, which 
can happen, for example, when rays reflect from a boundary 
or cutoff (see Fig.|4|i. In this case, each sheet of the Lagrange 
manifold will contribute to the phase of the eikonal solution. 
Each sheet can be associated with an eikonal solution; in the 
case of the reflection from the cutoff, one sheet is associated 
with the incoming wave, and the other with the reflected wave. 
In this case, the integral for 0(x) remains the same, and an ad- 
ditional phase shift is added to the outgoing wave. The value 
of this phase shift depends on the geometry of the Lagrange 
manifold through a quantity called the Keller-Maslov index, 
which will be discussed in more details below. 

The formula for the amplitude A (x) of the approximate so- 
lution is then written in terms of the determinant of the Ja- 
cobian matrix J(x) which arises when transforming from the 
coordinates (t,j3) to the physical coordinates x: 
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When the coordinate transformation becomes singular, the as- 
sociated field amplitude goes to infinity. This breakdown in 
the approximation for the amplitude is due to a singularity 
in projecting the Lagrange manifold from phase space to x- 
space. It is often the case, however, that this singularity can be 
"repaired" by finding a local solution by using the eikonal ap- 
proximation in A:-space. This eikonal solution is then Fourier 
transformed back to x-space, giving a local solution. The in- 
coming and outgoing eikonal waves can then be matched to 
this local solution in order to obtain an approximate solution 
without the amplitude singularity at the caustic. We will first 
illustrate how this procedure works by performing it for a sin- 
gle A^^ mode, which reduces to a one-dimensional (ID) prob- 
lem. We will then give the results for a full 2D problem, with 
a lower hybrid beam reflecting from the cutoff. 



A. The eikonal solution: jc-space 

Consider a single A^, mode, E{x;N^). The eikonal solution 
for this mode is straightforward, since it is a ID problem. The 
solution in x-space involves three steps: 



200 
100 


-100 
-200 



1.5 



FIG. 4: The ray (x(t),A:^(t)) for the ID eikonal solution. 

1. Trace the ray, which gives {x{T),kx(T)) 

2. Compute eikonal quantities e(T), 0(t), and A(t) 

3. Interpolate these quantities onto a grid in x to form the 
solution 

This process will give two eikonal waves, one which ap- 
proaches the caustic and one travels away from the caustic. 
These two waves will be matched to the local solution to give 
the complete eikonal solution. The two waves are identified in 
the algorithm by computing the sign of dx/dT along the ray; 
it is negative for the incoming part of the ray, and positive for 
the outgoing part of the ray. 

While this model is simple enough that explicit formulas 
can be written for the eikonal solution, we computed the solu- 
tion numerically since the point of this calculation is to com- 
pare numerical algorithms. In Fig. 2] the computed ray is 
shown. Fig.|5]shows the eikonal amplitude, phase, and polar- 
ization as computed along the ray. Notice that the amplitude 
blows up at the caustic near t = 6000. 



B. The local solution: eikonal in fci-space 

The eikonal solution in x-space diverges at the caustic, but 
this divergence can be "repaired" by estimating the solution in 
A;v-space. Here, the eikonal approximation can again be made: 



E(^,,f)=A(fe,)e(fe,)e'®'*' 



(25) 



This approximate solution is valid even when the ray is near 
the caustic point at kx = 0, since projecting from phase space 
to kf-space is well defined. The projection singularity encoun- 
tered in x-space is not a problem in this representation. 

First, consider the amplitude and polarization of this ap- 
proximation. These quantities can be evaluated along the ray 
as functions of T, and then interpolated onto k^ to get the over- 
all amplitude and polarization as functions of k^. The equation 
for the amplitude as a function of k^ takes the same form as 
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FIG. 5: The eikonal amplitude, phase, and polarization vector com- 
puted along the ray. Notice that the amplitude blows up at the caustic 
near T = 6000. In the plot of the polarization, the solid line is e.i (f), 
and the dashed line is £,{1:). 



Eqn. ( l23T l, with taking the place of x. This gives 



A{k^)e{k^) =Ao 



dkx 



-1/2 



e(^,(T),x(T)). (26) 



For the case A', = that we are considering, there are two 
components of e which are nonzero. This means that gener- 
ically, we need to construct two eikonal functions, one for 
each nonzero component of the polarization. However, we 
can use the wave equation [Eqn. and the dispersion ma- 
trix in Eqn. (|3]l to write a simple relationship between Eji{x) 
and E^{x): 

£x(x) ^ ^^^%2) ^^U^)- (27) 

So, we only need to construct the eikonal approximation for 
E^{x), and then take a derivative to obtain the approximation 
for Ex{x). Therefore, when considering the amplitude and po- 
larization A(A;x)e(fex), we are primarily interested only in the z 
covivpovitvii, A{kx)e,{k^) . Analytically, we expect this product 
to be constant, and when it is computed from the ray tracing 
data, it is constant to within a small numerical error 
We now have 



■'&(k. 
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where Aq is a complex, constant amplitude, and 



®{k,)= / x(4X. 



(29) 



As can be seen by the shape of the ray in Fig. |4] we can 
make a good approximation to the curve x{k:^) by assuming 
it is quadratic: 



x{kx) KiXQ + ykl, 



(30) 



where the constants xq and 7 are computed numerically, and 
found to be 



xo = 0.8747 and 7=3.137x10" 



(31) 



Some care must be taken with the sign when integrating to get 
@{kx), since the starting point of the integral is the starting 
point of the ray, which is a large positive value of k^. This 
gives an overall minus sign when doing the integral. Putting 
all of this together, we get, up to an overall complex constant 
amplitude, 



E-SK) = exp[i {-XQk^ - ^yk])] 



(32) 



Now, to get the solution in x space, take the Fourier trans- 
form of this function. The result can be written in terms of the 
Airy function, giving the local solution 



£,(x)=AoAi(-7-i/3(x-xo) 



(33) 



where Aq is a complex normalization that will be set by match- 
ing to the incoming wave. We can now use Eqn. ( |27] ) to find 
the X component of the local solution, and we get 



(34) 



C. The matched solution 

We now have numerical solutions for the incoming and out- 
going eikonal waves, which we can match to the local solu- 
tions in Eqns. (l33T l and ( l34l l to obtain the complete eikonal 
solution. In order to obtain the correct matching between the 
local and eikonal solutions, we will need to calculate the phase 
shift between the incoming and outgoing waves. This phase 
shift for our ID problem is given by the equation 



where pL is the Keller-Maslov index 



u = inerdex I 



inerdex 



dx 



(35) 



(36) 



and inerdex(M) is the number of negative eigenvalues of the 
matrix M (for a more complete discussion of the Keller- 
Maslov index in the context of constructing approximate 
fields, see ifToll ). For our problem, the slope of x{kj) is positive 
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for the incoming wave and negative for the outgoing wave, so 
we have ju = 0— 1 = — 1. 

It now remains to combine the incoming and outgoing 
waves in such a way as to form the complete eikonal solu- 
tion. This can be achieved by combining them with a switch- 
ing function, so that the local solution is used in the region of 
the caustic, and the eikonal solutions are used away from the 
caustic. The code developed for this paper uses a hyperbohc 
tangent function to achieve this switching: 

E(x) = [1 - w(.x)]Eiocal(-x) +H'(x)Eeikonal(-x), (37) 



where 



1 +tanh(^27"'/^(x-jc,„) 



(38) 



The matching point is taken to be a point where both the 
local and eikonal solutions are valid. The length scale over 
which the switching occurs is set using 7^'^^, since that is a 
natural length scale in the local solution. The complex ampli- 
tude of the local solution was set so that it matches the eikonal 
solution atx = Xm- 



■^z, eikonal 



Ai{-^y^{x„,-xo)) ' 



(39) 



The eikonal field E(jic) has now been constructed, and gives 
an approximate field which is good even in the region of the 
cutoff (see Fig.|6]l. 



D. The eikonal solution in two dimensions 

When computing the eikonal field in two dimensions, there 
are two important aspects which make the calculation more 
involved than the ID calculation. First, instead of simply trac- 
ing one ray, and using that to compute the eikonal field, a 
family of rays must be traced out (Fig. IT). Information from 
this family of rays is then used when computing the eikonal 
quantities 0(x), A(x), and e(x). This is for the most a part rel- 
atively simple change to implement numerically. Care must 
be taken, however, that the initial amplitude and phase at each 
ray is set appropriately. The initial phase is set by integrating 
(k • dx) from a chosen reference ray to the initial point of each 
ray in the family. The initial amplitude is somewhat more ar- 
bitrary. In this problem, we choose to set the initial amplitude 
for each ray using a gaussian profile, so that our resulting field 
will hopefully match the analytical solution shown in Fig.|3] 

The second critical distinction between the ID and 2D cases 
involves the question of how to find the local solution at the 
caustic. This is arguably the most difficult aspect of con- 
structing eikonal solutions in two dimensions. However, in 
our problem the simple geometry makes the calculation of the 
local field much easier. The x-dependance of the local solu- 
tion will be identical to that found for the ID case. The z- 
dependance is the simple plane wave exp{ik^z) multiplied by 
the gaussian amplitude profile which comes from the initial 
conditions (see local solution in Fig.fTOli. 
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FIG. 6: The eikonal and local solutions are combined using the 
switching function to obtain the matched solution. 



For this problem, the most computationally expensive part 
of the 2D calculation was interpolating the amplitude and 
phase data from functions specified along the rays, to func- 
tions of Figure [8] shows how the projection of the ray 
data onto the (x,z) -plane forms a nonuniform mesh. The 
eikonal data is specified on this mesh, and must be interpo- 
lated onto a uniform mesh in {x,z) in order to construct the 
eikonal solution. As an example. Fig. |9] shows the interpo- 
lated amplitude of the incoming eikonal wave. 

Once the eikonal phases and amplitudes of the incoming 
and outgoing waves are interpolated onto the (x,z)-plane, we 
are ready to match them to a local solution, and obtain our fi- 
nal fields. The same calculation as performed in the previous 



FIG. 7: (Color online) The family of rays for the 2D case, with pro- 
jections onto the {x,z)-, {Nx,z)-, and (.Y,A'v)-planes. 



FIG. 9: (Color online) Eikonal amplitude of the incoming wave, after 
interpolation onto (x,z)- 
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FIG. 8: Projection of the family of rays onto the (x,z)-plane. The 
eikonal data is calculated on this mesh. 



FIG. 10: (Color online) Absolute value of the z component of the lo- 
cal solution used for matching the eikonal waves in two dimensions. 



section gives the phase shift of k/2 for the outgoing wave. 
Matching the eikonal waves to the local solution (Fig. [TOl i 
which was described above, gives us our final results, shown 
in Fig.[TT] The result looks quite similar to our analytical so- 
lution; a detailed comparison will be made in SectionFVl 



IV. SEMICLASSICAL WAVE PACKET DYNAMICS 

In the previous section, the standard eikonal field approx- 
imation was computed for the lower hybrid cutoff model. 
While the eikonal approximation gives good results for this 
problem, it is somewhat unwieldy as a numerical algorithm, 
particularly because of the special care which must be taken 
in order to get a good field approximation at the cutoff. The 
process of finding a local solution, and matching it to incom- 
ing and outgoing waves, was fairly straightforward for our 
slab model, but more complicated geometries would make this 
process much more difficult. 

The standard approach of the previous section is not the 
only algorithm for constructing eikonal approximations, how- 



ever. For example the beam tracing approach described in 
Ref. |fl7|]. has been applied with some success to problems 
in plasma physics. fill In AMO physics 1 15l[T9i - l2lll and geo- 
physics [22] various methods using the idea of coherent states 
have been developed. The advantage of such methods is that 
they can naturally be interpreted in terms of the phase space 
dynamics of coherent states, which are a particular type of 
wave packet. Because the approximations being used for these 
methods are based in phase space, they have a natural advan- 
tage over configuration-space methods such as the standard 
eikonal approach described in the last section. Specifically, 
the presence of caustics is due to singular projections from 
phase space to configuration space, and since phase space 
methods do not require such projections, caustics do not arise 
in these methods. While the beam tracing method of Ref. fv^ 
can avoid some caustics, it is however still based on construct- 
ing a local coordinate frame in configuration space. This can 
lead to restrictions on the applicability of the method (e.g., 
the beam injection angle limit described in Ref. ifisll ) which 
would not be present for methods based in phase space. 
In this section, we will use a technique which was devel- 
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FIG. 1 1 : (Color online) Absolute value of Ex (top) and £7 (bottom) 
for the matched eikonal solution. 

oped as a semiclassical approximation for quantum problems, 
and which uses the semiclassical dynamics of wave packets 
for constructing approximate fields [151. This algorithm has 
a distinct advantage over the standard eikonal field approxi- 
mation, because it can automatically construct the field in the 
vicinity of a caustic, with no singularities in the solution. The 
amplitude of the solution remains finite at the cutoff, and the 
phase shift between the incoming and outgoing waves is com- 
puted automatically. This makes the wave packet approach 
ideally suited for implementing as an eikonal field approxi- 
mation algorithm. 

A. Description of the algorithm 

The algorithm described in this section for field reconstruc- 
tion was derived for solving quantum mechanical problems 
using the semiclassical dynamics of wave packets. As de- 
scribed by Littlejohn in Ref. |15], the dynamics of a wave 
packet can be approximated by considering ray trajectories in 
phase space. A wave packet initially centered at the point xq 
with the carrier wave vector ko will evolve into a wave packet 
centered at x{t) with wave vector k(f), where (x(f),k(f)) are 
the coordinates of the ray in phase space. The width and ori- 
entation of the gaussian wave packet then depend on nearby 



rays. If the nearby rays are diverging, then the wave packet 
will spread. Convergence of rays leads to focusing of the wave 
packet. It turns out that information about nearby rays can be 
approximated by a symplectic matrix S(f), and that the dy- 
namical equation for S(f ) has been derived. ifTsIl 

It will be convenient to group the phase space coordinates 
together as ^ = (x,k), with x = {x,z) and k = {kx,k-). With 
these definitions, there are three sets of ODEs that need to 
be solved in order to approximate the dynamics of the wave 
packet (for the derivation of these equations, see Ref. lilSIl ): 

= i(k.x-x.k)-D(<^(f)), (40) 

1 = J-VD(^(f)), (41) 
S = J-VVD(^(f))-S. (42) 

In these equations, is the eikonal phase, D is the dispersion 
function, and J is the symplectic matrix 

The gaussian wave packet that is following the ray ^{t)is then 
given by 

E(x,f) = ' \ > exp /k(f )-Ax(0 , (44) 

^det(A + iB) 

where the gaussian envelope of the packet, and the curvature 
of the phase fronts, are given by 

0(f) = ^Ax(f)-(D-/C)(A + /B)-^-Ax(f), (45) 
and where 

Ax(f) =x-x(f). (46) 

The matrices A,B,C and D are the submatrices of the sym- 
plectic matrix S: 

Note that x is the point where the field is being evaluated, 
while x(f ) is a point on the ray. The polarization e(f ) is the 
unit zero-eigenvector of the dispersion matrix evaluated at the 
point <^ (f ) on the ray, and Eq is the initial amplitude. Gener- 
ically, the matrix S (and thus its submatrices A,B,C and D) 
must be computed numerically. 

There are two important features of Eqn. (l44l) for the pur- 
pose of implementing it numerically. First, like other ray 
tracing algorithms, solving the wave equation has been re- 
duced to solving the ODEs in Eqns. (|40] | and (l42l i. [We solve 
Eqns. ( l40l i and (l4Tl i using standard integration algorithms, and 
Eqn. ( |42] | using the algorithm described in App.lAll The sec- 
ond important feature of Eqn. ( l44b is the one which makes 
this method very attractive for numerical implementation. It 
has been shown[23] that the matrix A + /B is never singular. 
This means that its determinant is never zero, and thus the 
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wave packet amplitude never diverges. This is true even at the 
location of caustics, where ordinary eikonal methods would 
break down. Additionally, the square root •\/det(A + /B) is a 
complex function. By choosing its phase to be continuous in t, 
the phase shift that the wave undergoes at caustics is automat- 
ically calculated. Because Eqn. (l44l l is vaUd and well behaved 
at caustics, it is ideally suited for implementation as an numer- 
ical algorithm for wave field calculations. Such an algorithm 
should be able to correctly reconstruct the wave field, even 
near caustics. Because of this, it would be superior to the stan- 
dard eikonal method described in the previous section, where 
the caustic at the cutoff required special handling. 

Fourier transforming the wave packet in Eqn. (l44b from 
time to frequency gives a time-independent wave field: 




200 t 400 



Eco(x) = / e"^'E{x,t)dt. 



(48) 



While there is no reason in general that this integral should 
converge ||24|, integration over a finite interval in time gives 
a field which can be interpreted as a "beam" of finite length. 
In the next section, we show the results of implementing this 
approximation for the lower hybrid wave reflecting from the 
cutoff in a plasma. 



B. Setting initial conditions 



FIG. 12: The initial value of for the wave packet was set so that 
the packet hits the cutoff "head-on," i.e., = (solid line) at the 
same time that dx/dt = (dashed line). 



and its relation to realistic boundary conditions, is an area for 
future research. 

Numerically, a parameterization such as that described in 
Ref. 126] can used to set Sq. However, after exploring various 
initial conditions, it was found that a relatively simple param- 
eterization could reproduce our analytical solution fairly well: 



Perhaps the most difficult issue to address when using the 
wave packet approximation described in Sec. lIV Al is the ques- 
tion of how to set the initial conditions for the wave packet. 
The initial point of the ray sets the initial location of the peak 
of the wave packet, but the packet shape is set by the initial 
value of the matrix S. Fortunately, this initial value can be 
set even after S(f) is calculated, because of the form of the 
solution. Eqn. (l42b implies that a solution S(f ) has the form 

S(f)=exp{J.VVD[^(f)]}-So, (49) 

where So = S(0) is the initial condition for the matrix S ll25ll . 
Since the solution has this form, another solution with differ- 
ent initial conditions So can be constructed simply by multi- 
plying on the right: 

S(f) = S(f)-Soi-So. (50) 

For our calculations, we used the identity matrix as the initial 
condition; So = id. Then, after finding S(f), we applied new 
initial conditions which gave the desired packet shape. 

The symplectic matrix So of initial conditions is a 2n- 
dimensional matrix, where n is the number of spatial di- 
mensions. The space of symplectic matrices has dimension 
«(2« + 1) = 10. For a generic problem, the constraints used to 
set So would be related to the source of the wave, such as an 
antenna or a waveguide. In the problem considered in this pa- 
per however, we are studying the reflection of a beam from the 
cutoff and so the details of the source of the beam are not as 
important. Therefore, we are primarily interested in choosing 
So so that the resulting beam has a shape which matches our 
analytical solution. The choice of So for a generic problem. 



( (gV' )' ^^^^ 

where G is an invertible n xn matrix. For our problem, we 
found that a real-valued, diagonal matrix form for G gave us 
reasonably good results: 

The parameter (7, was estimated by taking the analytical so- 
lution shown in Fig. [3] and fitting a gaussian function of z to 
the incoming beam. This process is similar to that which was 
used to set the initial amplitude for the eikonal solution in the 
previous section. Setting the parameter Gx is somewhat more 
difficult. Numerically, it was found that varying could sig- 
nificantly change the spreading of the wave packet. Fairly 
reasonable results were found by setting Gx in the following 
way. 

Consider the gaussian envelope of the wave packet at time 
f . Contours of the wave packet are ellipsoids which make an 
angle 0(f) with the z-axis. The value of Gx was chosen so that 
the wave packet hits the cutoff "head-on." This was done by 
numerically finding the value of Gx such that 0{t*) = 0, where 
f * is the time when the ray reflects from the cutoff. This value 
can be found visually by plotting 9{t) and dx/dt together. At 
the cutoff, dx/dt = since the ray reflects in x. So, we varied 
Gx until 9{t)—Q at the samepoint that t/.ic/iif =0 (see Fig.[T2]i. 
Numerically, this results in the values Gx = 0.1174 and G^ = 
0.2301. 
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FIG. 13: (Color online) Ellipses drawn at the l-CT amplitude level 
show the time dynamics of the wave packet. Notice the reflection 
from the cutoff, and the spreading of the wave packet. 
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FIG. 14: Wave packet amplitude, A = |det(A + (B)|^'''^ asafunction 
of ray parameter t. The amplitude increases at the cutoff and then 
decreases as the wave packet spreads. 



C. The wave packet solution 

Having set the initial conditions as described above, the 
wave packet solution was computed, and the fields were cal- 
culated using Eqn. (l48T l. Figure [13] shows l-a contours of the 
wave packet's gaussian envelope as a function of time. Fig. [14] 
shows the wave packet amplitude, A = | det(A + /B) | ^ 'Z^. No- 
tice how the amplitude gets large at the cutoff, but does not go 
to infinity as the eikonal amplitude does. 
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FIG. 15: (Color online) Locations in the (x.z) -plane where the com- 
parisons are performed are marked by white lines. The E~ component 
of the field is compared atx = 0.905, and and at z = 0. 



Rather than comparing the entire 2D fields, we compare the 
fields at fixed x (and z) position. Figure[T5]shows the locations 
where these comparisons are taken, while the actual compar- 
isons are shown in Figs. [16] and [17] As can be seen in the 
comparisons, the error in the eikonal construction (calculated 
by subtracting the approximate field from the analytical so- 
lution) is about 3%, while the error in the fields constructed 
using the wave packet is about 10%. Looking at the shape of 
the error in Fig. [17] it is not unreasonable to suspect that the 
spreading of the wave packet is causing the error to be this 
large, since the errors are largest at the edges of the beam. 

The spreading seen in the semiclassical wave packet dy- 
namics would be an interesting effect to examine in future 
research. The wave packet dynamics uses second derivatives 
of the dispersion function, and so in a sense incorporates more 
information about the physics of the problem than the standard 
eikonal method. It would be interesting to study how much 
of the spreading is due to numerical effects, and how much is 
real physics. Another aspect of the wave packet method which 
would be interesting to study would be the setting of the ini- 
tial condition So for the wave packet. Numerical exploration 
of various values of So showed that the spreading can change 
quite a lot depending on the choice of So. Perhaps a differ- 
ent choice of So could even reduce the errors seen in Figs. [16] 
and [it] Further study of this aspect of the problem would be 
needed to determine if this were possible. 



VI. CONCLUSION 



V. COMPARISON OF SOLUTIONS 

In Sec.|ll] the analytical solution was found for a lower hy- 
brid wave reflecting from a cutoff in a slab model. The stan- 
dard eikonal field approximation was constructed in Sec. [Ill] 
and in Sec. lIVI an alternative algorithm based on the semiclas- 
sical dynamics of wave packets was used to obtain an approx- 
imate solution. In this section, we compare the two approxi- 
mate solutions to the analytical solution, in order to estimate 
their accuracy. Note that the eikonal solutions were normal- 
ized to match the analytical solution. 



In this paper the reflection of a lower hybrid wave from a 
cutoff was studied, in order to illustrate two different eikonal 
methods for constructing approximate solutions. Comparing 
these methods in this simple model allowed us to examine the 
relative merits of each. 

The standard eikonal techniques used in Sec. [Ill] give very 
good results, but special care is needed near the cutoff so that 
the field near the caustic is calculated correctly. This is rela- 
tively straightforward to do in the slab model discussed here, 
but in more realistic geometry it would be much more diffi- 
cult. Also, the calculations needed for calculating the fields 
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FIG. 16: Comparison of approximate solutions to the analytical so- 
lution. This is a comparison of the real part of computed at 
X = 0.905, showing a slice across the beam. The top panel shows 
the exact solution, the middle (bottom) panel shows the error in the 
eilconal (wave packet) field approximation. 

near the caustic are difficult to automate on a computer, and 
thus limit the usefulness of this method for generic eikonal 
calculations. 

The algorithm reported in Sec. |IV]is based on the semi- 
classical dynamics of wave packets. While this algorithm was 
developed for quantum mechanical calculations, it is also ap- 
plicable to wave problems in plasma physics. Because the ap- 
proximations used in this method are based on consideration 
of nearby rays in phase space, it does not have the projection 
singularities at caustics that the standard eikonal method has. 
This means that the fields based on this method do not blow up 
at the cutoff. More importantly from a numerical viewpoint, it 
means that no special treatment is needed in order to compute 
the correct fields at caustics. Because of this, one would ex- 
pect that this algorithm would be much better suited for com- 
puting eikonal field approximations than the standard method. 
There are a few downsides to using this method, however, 
which would need to be addressed in order to create a robust 
code based on this algorithm. In particular, the issue of how to 
initialize the wave packet in order to get a good match to ini- 
tial or boundary conditions would have to be addressed. The 
method used in this paper turned out to be sufficient for the 
simple problem considered here, but a better method would 



FIG. 17: Comparison of approximate solutions to the analytical so- 
lution. This is a comparison of the real part of E,, computed at z = 0. 
The top panel shows the exact solution, the middle (bottom) panel 
shows the error in the eikonal (wave packet) field approximation. 



need to be developed for more generic problems. There is 
also the question of what to do when the spreading of the 
wave packet becomes too large, causing the nearby ray ap- 
proximation to break down. This would be an issue especially 
in problems where the system size is relatively small, since 
smaller wave packets tend to spread faster. It could also be 
an issue in problems where very long ray trajectories would 
need to be followed. Despite these shortcomings, this wave 
packet algorithm holds promise for performing computations 
in a variety of realistic plasma physics problems. 

The calculations reported here show that it is possible to 
develop a code to construct approximate wave fields based 
on the semiclassical dynamics of wave packets. With care 
taken to incorporate the effects of antenna geometry, better 
plasma models (including damping), and more reaUstic ge- 
ometry (perhaps even three dimensions), such a code could 
produce useful and interesting results. 
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Appendix A: Solving for S(f) 

The symplectic matrix S(f) describes the dynamics of 
nearby rays, and gives the parameters of the gaussian wave 
packet used for constructing an approximate field in Sec.lIVI 
The symplectic nature of this matrix is an important property, 
and it is not preserved by a generic ODE solver. We used 
the following algorithm for solving Eqn. (l42l) . which has the 
property that S(f) can be kept symplectic to desired order in 
At. 

Rather than numerically solving Eqn. ( l42b for the compo- 
nents of S(f), first notice that the solution has an exponential 
form: 

S(f + Af ) = exp Q^'^^' J • VVD{z{t'j)dt'^ • S(f) (Al) 
~exp(J- VVD(z(f))Af)-S(f). (A2) 



This approximation for S after a time step At has a par- 
ticularly nice property: it can be shown that the matrix 
exp(J • VVD(z(f))Af) is symplectic. The product of sym- 
plectic matrices is also symplectic, so this approximation for 
S(f + Af) is also symplectic. 

Now, in order to compute the exponential of a matrix, ex- 
pand the exponential above to desired order in Af: 

S(f + Af)^ (^id + AS + ^(AS)2 + ...^ .S(f), (A3) 

where 

AS=J- VVD(z(f))A/. (A4) 

While the expansion in Eqn. (IA3l l does not maintain the sym- 
plectic property which we desired, higher order terms can be 
included to yield an approximation which is closer to being 
symplectic. The error should go like (Af)"+', where n is the 
highest order kept in the expansion. For the results reported 
in this paper, n = 4. 
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